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1. Introduction. 

In this paper we consider a one-dimensional seismic inverse problem: 

Our goal is to construct a parameter estimation scheme for the hyperbolic 
model equations treated in [7] extending these ideas to allow the estimation 
of discontinuous coefficients (including the location of discontinuities). 

The approach taken here contrasts with that taken in [7] in that we consider 
a different decomposition of the wave equation, yielding a different operator 
and state space for theoretical arguments. We combine these ideas with a 
variation of the approximation scheme in [15] which was developed to treat the 
problem of estimating discontinuous coefficients in parabolic systems. 

The underlying theoretical approach to the identification problem follows 
the outline of other related papers (eg. [3], [4], [5], [6], [7], [15] in that we 
define the (infinite dimensional) identification problem and construct associated 
approximate identification problems; under compactness assumptions on our parameter 
space, we prove convergence of parameter estimates and approximating state 
variables. Our arguments here are based on an application of the Trotter-Kato 
Theorem. In Section 5 we describe the numerical algorithm and conclude with 


some preliminary computational examples. 

We shall employ standard notation throughout, using, for example, H p (n) 

and Wp(n) to denote the usual Sobolev spaces on n. If n is not specified, it 

is assumed that n = [0,1]. Further, given any w e L°°(0,1) satisfying 

0 < w < w(x) < w, for almost all x, we shall define the -j- - weighted H^(0,1) 

-i w 

0 f 1 

SDace, denoted H (w), with inner product <u,v> e 4_ uv ( anc | associated 

w j Q w 

norm |*| ). Similarly, given any c>0 we define the c - weighted real line by 
R(c) with inner product of two elements u, v e K(c) defined by cuv. Finally, 
throughout we shall write I : Y -*■ Y to denote the identity operator, where the 
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space Y will be clear from the context. 


2. The Identification Problem. 

As in [7], we consider the problem of estimating unknown parameters 
that appear in the system 


A 


( 2 . 1 ) 


— (E(x) |^-) , t e (0,t] , x e [0,1] , 

at p (x) ax 3X 


— (t,0) - k^Ct.O) = s(t;k) 


at 


(t,l) + k 2 |^(t,l) = 0 


\u(0,x) = 4> (x) , f£(0,x) = ip(x) 


where, for the seismic problem, u represents subsurface particle displacement 
resulting from a disturbance at the earth's surface (x = 0). We model the 

Or "I 

disturbance here using the source term s(t;k) (s e H (0 , t ) ) which may involve 
an unknown vector parameter k. Other parameters of interest include p and E 
(p > 0, E>0), which represent the density and elastic modulus, respectively, 
of the earth, and k-j and k 2 (k-j > 0, k 2 >0), which occur in the elastic boundary 
condition at x = 0 and the absorbing boundary condition at x=l, respectively. 

A more thorough discussion of the use of model equations (2,1) in the context 
of the seismic problem may be found in [7]; similar models are treated in [9] 
or [11], where boundary conditions are imposed instead on a semi-infinite strip. 

The estimation problem of interest then is to identify the unknown 
parameters using observations of the system, which may be given in the form of 
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displacement, velocity, or pressure (stress) measurements (corresponding to 
u, u t , or Eu x in (2.1)). Throughout we shall assume that unknown parameters 
are expressed in the form 

9 = (q^*)* ( * ) » q 3 » ^5) = ( p ( . ) * £(*)» , k 2 , k) 
where q belongs to the admissible set Q, 


Q = f (q-, ,q 2 ,q 3 ,q4,q 5 ) e H° x H° x K 2+k 


for i = 1, 2: 0 < q . < q^x) <_ for almost all x e [0,1] 

and q.j is continuous at x= 0, 1 ; 0 < R3 £ q 3 ± q 3 l > 0} 


Due to the fact that the spatial domain [0,1] represents an inhomogeneous, 
layered medium, it is likely that the two spatially varying parameters, q-j 
and q 2 , will be discontinuous with locations of discontinuities corresponding 
to abrupt changes in subsurface structure. In keeping with the objectives 
of the seismic inverse problem, i.e., to assemble as much information as 
possible about the nature of the subsurface medium, our goal will be to 
determine the location of discontinuities, as well as to identify the spatial 
variation in q-| and q 2< 

Before considering a precise statement of the parameter estimation 
problem, we first turn to an abstract formulation of (2.1), which is necessary 
for the convergence theory we develop in Section 4. To this end, we (formally) 
rewrite (2.1) as a first-order system in v(t) = (v^(t), v 2 (t), Vg(t)) ^ 

(u(t,0), u t (t,*). q 2 (*)u x (t,*)) e F x H 1 x H 1 , namely. 
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0 I Q 0 

v(t) = | 0 0 q-,0 J v(t) 

0 qgD 0 


( 2 . 2 ) 


Vo(t) 


Q " q 3 q 2^°^ v l = 


(q 2 0)v 2 (t) + c * 4 v 3 (t) ) | ^ 


= 0 , 


v(0) - ( <J>(0) , <?(•)» q 2 (* ) 4> x (* ))T , 


where I Q : C[0,1] ■* ]R is defined by I Q w = w(0) and D represents the spatial 
differentiation operator. Such a decomposition of state variables can be found 
elsewhere in the literature (see for example [10] in a different context): It is 

a natural one for the application we have in mind because the components 
u(t,0), u^, and q 2 u x are quantities needed to express the total energy of this 
system. In addition, these quantities are often readily observable in 
practice and thus their inclusion as components of the state variable facilitate 
our study of the estimation problem. 

We use the form of (2.2) to define a parameter - dependent operator 
A = A(q) and associated Hilbert space X = X(q), where dom AcX and X is chosen 
in such a way to ensure that A is dissipative in X. In particular, for any 
value of q e Q we define X(q) = JUq^q^O)) x H^(q-j)x H^qg) (see Section 1 
concerning notation) with associated norm denoted by ||*|| . It is easily seen 
that a positive constant u = y(SLj» q^ ’> i = 1,2,3) may be found so that 
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(2-3) 1|UII « llz |] q < u||z |1 

for all qeQ, where ||»|| denotes the Rx H°x H° norm; it thus follows that the 
X(q) norms are uniformly equivalent as q ranges over Q. Due to the set-wise 
equivalence of X(q) for all qeQ, we can consider an element z in X(q) to be 
an element in X(q), for q e Q, or even consider z in R* H^x H®, and do so 
without a change in notation for z. 

The abstract differential equation associated with (2.2) may be stated 
precisely as 


(2.4) 


z(t) = A(q)z(t) + F ( t ; q ) , t e (0,t] , 


[z(0) = z Q (q) 


where the transformation (z-|,z 2 ,z 3 ) = (v.j,v 2 ,(v 3 + (x-l)s)) has been made in 
order to obtain homogeneous boundary conditions. In (2,4), z(t) e X(q), and 


A(q) = 


/ ° 'o 0 \ 


V 


0 q-jD 
q 2 D 0 


/ 


with dom A(q) = (z e X(q) fl (1R x x ) 


z 3 (0) - q 2 (0)q 3 z ] = 0 


q 2 (l)z 2 (l) + q^z 3 (l) = 0); the nonhomogeneous term (due to the transformation 
to zero boundary conditions) and initial condition are given by 
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1 

f 0 ) 


{ <0(0) \ 

F ( t ; q ) = 

-q-| (• )s(t;q 5 ) 

^ (*-l)s(t;q 5 ) j 

• z 0 (q) = 

♦ (•) 

yq 2 (*)<O x (*) + (- -1 )s(0;q 5 )/ 


respectively. 

Theorem 2.1 . A(q) is the infinitesimal generator of a Cq - semigroup of 
contractions T(t;q) on X(q); thus for each q e 0, there exists a unique mild 
solution of (2.4) with the representation 

ft 

(2.5) z(t;q) = T(t;q)z Q (q) + T(t-s;q)F(s;q)ds . 

The proof of the theorem follows from standard results from the theory of 
semigroups (see, eg., [17; pp. 14, 106]) once it has been verified that A(q) 
is densely defined and dissipative in X (q ) and that for some x>0, R(x I - A(q ) ) = X(q). 
It is easy to see that dom A(q) = X(q) since X(q) is equivalent to IRx * H^, 
while dissipativity of A(q) is a consequence of the topology chosen for X(q). 

Moreover, using standard results from the theory of two-point boundary problems, 
one can show that R(xl-A(q)) = X(q) for any choice of x>0. 

Remark 2,1 . Strong solutions of (2.4) (i.e., solutions where z Q (q) e dom A(q)) 
will satisfy (2.4) as well as (2.5). In this case, (v-|,v 2 ,v 3 ) = (z^ ,z 2 ,z 3 + (l-*)s) 
satisfies the formal equations (2.2) so that the identification may be made 
between z-j and u(t,0), z 2 and u t (t,*)» and z 3 + (l-.)s and q 2 (* )u x (t,* ) . 
where u is a solution of the original system (2.1). Hence, we shall develop 
a convergence theory based on formulations (2.4) , (2.5), while keeping these 
relationships in mind throughout. 
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Remark 2.2 . That our approach is advantageous from both a theoretical and 
computational viewpoint is now apparent: Physical principles demand that 

q 2 u x be continuous even though q 2 , and thus u x , are discontinuous. By defining 
our operator as above, we satisfy this continuity condition simply by requiring 
that the third component of the state variable belong to h"*(0,1). 


We turn now to a precise statement of the parameter estimation problem. 

The unknown parameter, q, of interest contains for its first two components the 
piecewise continuous coefficients q-j and q 2 ; we further parametrize q-| and q 2 by 


v 

q-|(x) = cx 0 (x) + J a.(x)H 5 (x) , 

( 2 . 6 ) 

V 

q 2 (x) = 6 0 (x) + j^B..(x)H c (x) , 


where H denotes the usual Heaviside function (H_ (x) = 1, x e (c,l], and 
(x) = 0 otherwise) on [0,1]. Using this representation, the problem of 
estimating q-| and q 2 is equivalent to the problem of identifying 2(v + l) 
continuous "pieces", <xq(x), ..., a v (x), Bq(x), ..., B v (x), and the location 
of v discontinuities, 5-j, ..., e . Given v, our goal then is to estimate 
P = (c*q> ® i » •••» » Bq» * •••» B^ , 5 -j , •••, q^» q 4’ ^5^ where p 

belongs to a prescribed constraint set PcP, 

f\j 

(2*7) P = { (otQ 9 • • * ) Ct^9 3q> •••* 3^5 £*|» •••» 9 Qg) 

6 ( 2< fj ) x R v+2+k 

V V 

(“o + ^ a i H r. » e 0 + J-j B i . » q 3’ q 4’ q 5^ 6 ^ * 
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Although our estimation theory will be developed for the fully parametrized 
set P, it is clear that any choice of p e P generates a physically meaningful 
entity q = q(.p) e Q through the construction (2.6) for q-j and qg, and that 
sequential convergence p J -*• p in the topology for P implies the associated 
parameters q J = q J (p J ), q = q(p) e Q satisfy q J -*■ q in the topology on Q. We 
shall often make use of this association in the sections to follow. 

The parameter identification problem may now be stated as follows: 

(ID) Minimize J(p) = \ ||Cz(t J .;q(p)) - y^ || 2 

J 

over p e P, where q = q(p) e Q, z(t;q) is the solution to (2.4), 

C: X(q) -*• F * is a (continuous) observation operator and 

y . e R x H° * H° denotes observed data corresponding to time t.. 

J J 

Remark 2.3 . The estimation problem defined above allows for observations of 
particle displacement (at the surface), velocity, or pressure. Although 
it is not explicitly stated, our theory also includes the case of spatially 
distributed displacement observations below the earth's surface; in this case, 

the observed quantity is given by z(t.;q) =<}>(•) + J z~(s;q)ds, and 

J Jo c 

observations are given by y., where y,-* z(t-;q) e H^. 

J J 0 

Remark 2.4 . The parameter estimation problem described above is an example 
of a large class of problems ("inverse" problems) that are widely known to be 
ill-posed (see, for example [13], [16])from both a theoretical and computational 
standpoint. We shall not address here the nature of the difficulties that 
arise, among them a lack of continuous dependence of parameters p on observations 
y. , and the general nonuniqueness of an optimal parameter p. However, it is 

J 
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appropriate to note that one may circumvent some of these inherent problems by 
guaranteeing a type of "problem stability" (see [2], [4], [11], and 
[12; Remark 5.1]). There are a number of ways in which such stability may be 
obtained for the parameter estimation problem: These include such techniques 

as regularization [13], embedding methods [1], or parameter set compactness 
assumptions [2], [5], As will be evident in later sections, we shall take the 
latter approach through the assumption that P is compact. 


3. Approximation. 

In this section we construct a spline-based framework for the approximation 
of the state variable z(t), where z satisfies (2.5). As these approximation 
spaces will be parameter-dependent, we shall assume that p e P is given, 
where, p = (a Q , a-|, Bq, B-j , £ , q 3 , q 4 , q g ) and 5 = £ -| e (0,1) (for simplicity 
and without loss of generality, we shall take v = 1 throughout). Using 
(2.6) to construct q^ , q 2 , we shall as usual associate with p the corresponding 
q = q(p) *e Q. 

N 

The construction of approximation spaces X (q) is as follows: For 

N 

k = 0, ...» 2N we define uniform spatial mesh points by s^ = k/2N and denote by 

II 

the k in (standard) basis element for the space of continuous piecewise-1 inear 
B-splines with knots at { s ^1 ; i.e., is characterized by SjJsj) = 6.^, 
j,k = 0, .... 2N,[18]. We then transform to a parameter-dependent basis through 
the invertible mapping g: [0,1] [0,1], 


(3.1) g(x) = 


I x/2? , 


S' 


[(X + 1 - 

and define bJJ(x) = sj!(g(x)). 


0 <_ x < 5 

25 )/2(l -? )» 5 1 x 1 1, 

'v w 

Thus sd{ B^J is the space of linear B-splines 
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with knots at xj^ e g _1 (sj^) (x£! = k£/N, k= 0, .... N; xj^ = 5 + (k- N)(l -5 )/N, 
k=N + l, 2N). The computational advantages of usings - dependent elements 
will be discussed in Section 5. Finally, we define approximation spaces 
X^(q) = sp { B^} where the basis elements B^, j = 0, 4N, are defined as follows 

(we use 0 to represent the identically zero function on [0,1]): 


Bj = (0, bJ, 0), j = 0, .... 2N-1; 

B 2N = (0 ’ • “ q 2 (1) ®2N }i 

Bj = (0, 0, B 4N _ d ) 9 j = 2N + 1, 4N-1; 

It is clear from the above construction that X^(q) <= dom A(q). It is also 
clear that if q t q then X N (q ) f X N (c} ) and, in fact, X f,| (q ) £ dom A(q ). 
Thus, as we iterate on q, we must take some care with the changing domains of 
the operators; we shall address this difficulty in the next section. 

We define approximating equations associated with (2.4) by defining 
operators A^(q) to "approximate" A(q). Here we take A *(q) = P^(q)A(q)P^(q) 
where P"(q) : X(q) -> X^(q) denotes the orthogonal projection (with respect 
to the topology on X (q ) along (X^(q) )> The differential equation on X^(q), 
equivalent to a system of ordinary differential equations, is given by 


(3.2) 


z N (t) = A N (q)z N (q) + P N (q)F(t;q), t e (0,t], 
^z N (0) = P N (q)z 0 (q) , 


where z^*(t) e X^(q). Using standard arguments (see, e.g., [7]) it is easy to 
show that A N (q) is bounded and generates a Cg - semigroup of contractions T N (t,q) 
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on X(q) ; thus, there exists a unique mild solution z N (*;q) e C(0,t;X^(q)) of 
(3.2), expressed as 


.N/ 


rN/ 


,N/ 


ft 


rN/ 


,N/ 


(3.3) z (t;q) = T ( t ; q ) P (q)z Q (q) + T (t-s;q)P (q)F(.s;q)ds 


for t e [0,t], 

The approximate identification problem associated with (3.3) is the 
following: 


(ID N ) Minimize J N (p) = \ || Cz N (tj*,q(p)) - y^l 2 

J 

N 

over p e P, subject to z ( * ;q (p ) ) satisfying (3.3). 

It is not difficult to argue, using the matrix representation for (3.2) 

(see Section 5), that, for z e ]R * * H°, the mappings p-> P N (q(p))z and 

p T M (t;q(p))z are continuous, the latter uniformly in t e [0,t]. Under 
reasonable assumptions of continuity on the mappings q 5 + s(t;q 5 ) and 
qg -> s(t;qg) (again uniform in t) we also obtain the continuity of p -»• ZQ(q(p) ) , 
p -*- F(t;q(p)). We thus obtain, from the continuity of p J (p), the following. 


Theorem 3.1 . Assume P is compact. Then for each N there exists a solution 
p N e P to ( ID N ) . 

4. Convergence. 

In this section, we will establish that a subsequence of the parameter 
estimates generated by solving the approximating identification problems (ID^) 
does indeed converge to a solution for the original identification problem. 

(We note that in practice, we and others working with similar schemes, have 
observed direct convergence of the estimates; in fact, under the assumption 
of a unique solution of (ID), full sequential convergence is guaranteed.) 

We will use the theoretical framework developed in [5], and discussed 
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in the context of the seismic problem in [7]; we refer the interested reader 
to these references for proofs and details, and here simply state relevant 
results without proof, focusing instead on the details which are new to our 
formulation of the problem. 


Throughout this section we will assume P is compact and that an arbitrary 

N M r\j 

sequence {p } in P has been given with p p e P, where 



and £ satisfies s < £ < 1 - <5 for some 6 > 0. We recall that the associated 
parameters = q^(p^) and q = q(p) in Q then satisfy q" -*■ q in Q (i.e., in 


H°x 


H° x P 2+k ). 


An important intermediate result in our convergence theory is to show that 
z^(q^) z(q) in X(q), where z^(q^) is a solution of (3.3) and z(q) is a 
solution of (2.5). To this end, we first demonstrate convergence, in an 
appropriate sense, of the semigroups. We will use the following version of 
the Trotter-Kato Theorem (see [5]). 


Theorem 4.1 . Let (B,|*|) and (B^,|*|„) f N= 1, 2, ..., be Banach spaces and 

let : B B^ be bounded linear operators. Further assume that T(t) and 

N N oj 

T (t) are Cq - semigroups on B and B with infinitesimal generators A and A , 

respectively. If 

(i) lim |it N z| n = ) z | for all z e B , 

N 

(ii) there exist constants M,w independent of N such that 

|T N (t)| N _< Me l ° t , for t >_ 0, and 

(iii) there exists a set V c B, P<=dom A, with (Xq-A)P = B 
for some Aq > 0, such that for all z e V we have 
|A^ir^z - TT^zj^ + 0 as N ->■ 


oo 
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then | T ISf ( t )-rr N z - ir N T(t)z| N + 0 as N + », for all z e B, uniformly in t on compact 
intervals in [0,«). 


We shall use the above theorem with B = X(q) and B^ = X(q^). For N=l, 2, ..., 
the linear operators ir N : X(q) -> X(q N ) are defined by ir N (z 1 »z 2 »z 3 ) = (tt^z-j .T r^.i^Zg) , 
where the positive constants irlj 1 are defined by ^3) » 

it 2 = qJ^O)/^ q^^» and 11 3 = An essential property of n N (which can 
be easily verified) is that maps elements of dom A(q) to dom A(q^). It is 
clear that 1 as N -»• ~ for i = 1, 2, 3, so that for z e X(q), 

(4.1) ||U N -I)z || N ->• 0 as N + » , 

and that Theorem 4 . 1 ( i ) is satisfied; further, the operators are bounded uniformly 
in N. We take T(t) and T^(t) to be the Cq - semigroups generated by A(q), A N (q N ), 
respectively; condition (ii) of the theorem follows from the fact that for all 
N, is a contraction semigroup. As a first step in the verification of (iii), 
we select V = dom A(q)fl (R * W^ * w^) and again appeal to the theory of two-point 
boundary value problems to claim that, for any X > 0, (xl - A(q) ) maps V onto 
R x L°° x L", which is dense in X(q). To satisfy the hypotheses of Theorem 4.1 it 
only remains to show that for z e P, we have ||A^(q^)ir^z -n^A(q)z|| 0 

q r 

as N “. 

It is helpful to establish several lemmas, the first of which is a 
generalization of standard estimates regarding the usual linear spline 
interpolating operator, i^, defined on the mesh {xj^, k = 0, ...» 2N} (i.e., 

it 2 N M r\j » 1 

i t = T f(xHB^, in the notation of Section 3). The corresponding state 
k=0 K K 

space interpolation operator of interest is the mapping 

I^:(IRx x ) -> (IR x x ), where we define lx i^ x i*'*. 
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N 

Lemma 4.1. There is a constant Cg independent of N and q such that for any 
z e dom A(q), 

(4.2) || I N z - z |[ qN < CgN -1 1 | D(I N z - z)|| qN , 
where we define DCz-j^.Zg) = (0,Dz 2 ,Dz 3 ). Further, 

(4.3) || D(I N z - z) || M + 0 as N + •. 

q IN 

Proof: For any f e H 1 , it is an established result ([19], eqn. (2.16)) that 

(4.4) |i N f-f| 0 l ( ir N ) ~ 1 1 D(i^f - f ) | gj 

in addition, one obtains | D(i N f - f)|g->- 0 as N •+•« from [19], eqns. (2.10) and (2.17), and 
the dense inclusion of H in H . Both (4.2) and (4.3) follow immediately 
from these estimates and the statement of norm equivalence (2.3). |~| 

N 

Lemma 4.2. There exists a constant c^ independent of q and N such that, 
for any z e X N (q N ), 

(4.5) || Dz || q N £ c-j N|| z || N , 
for N sufficiently large. 

Proof: From (2.3), || Dz || 2 N < y 2 | Dz 2 | 2 + y 2 | Dz 3 | 2 where, for j = 2,3, 
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2 


dx < 12 


2N N 
I (X - 
k=l K 




•x 

N 

k- 



2 


dx , 


from the Schmidt inequality [19; p. 7]. Since s’ satisfies 6 < s' < 1-6, it 

N N 

is clear that, for N sufficiently Targe, £ satisfies 6 <_s <_ 1 -6 and that 
x k - x k-i 1 «/N. Thus, ||Dz ||* N < 12(uN/<$) 2 (|z 2 | 2 + |z 3 | 2 ) and, using (2.3) 

again, the result in (4.5) is obtained. |~| 


Lemma 4.2 is a generalization of the Schmidt inequality for our product 
space and topology. The next lemma is such a generalization of the First 
Integral Relation for linear splines (i.e., |D(i N f)| Q £ | Df | Q for f e H 1 [19; p. 16]) 
the proof in our case follows directly from that relation and the statement of 
equivalence of norms (2.3). 

N 

Lemma 4.3. There is a constant c 2 independent of N and q such that 
(4.6) ||DCr N z)|| N iC 2 ||Dz|| N 
for any z e X(q), and any q e Q. 

The spline estimates in the preceding lemmas may now be used to obtain 
convergence of P^(q^) to I, in an appropriate sense. To simplify notation 
we shall write e P^(q^), where no confusion exists. 



16 


.'v. 


Lemma 4.4 . For z e dom A(q) 


(4.7) ||(P N -I)„ N z|| „ < || T N .„ N z - , N Z 


<-|(||(, N - I)Dz|| 2 n +||D(I N z-z)|| 2 n )’ s 

q q 


N 


for a constant c 3 independent of N, q and z; moreover, for z e X(q), 


(4.8) || (P N - I )7r N z|| + 0 as fU « . 


Proof: For any z e dom A(q), the definition of ir N guarantees tt N z e dom A(q N ) 


and hence I^ir^z e X^*(q^). Thus 


ll(p N -i)A|| q 2 N < ihi n -i)A|| 2 n 


2 / | . N / N % / N % .2 

<_y (|i U z ) 2 " z ) 2 1 0 


4. |*N/ N \ , N \ ,2x 

+ |l (n Z) 3 - (tt Z) 3 |q) 


From (4.4), it follows that 


i N (A)j - £ («N)-‘ |D(iV'z)j - U"*>j)IS 


N_\ ,2 


-2 , N.x m2 


for j = 1,2, where 
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|D(1 H (» N z)j-(» N ;Oj)l 0 £ |D(1 , Vz) j -1 , ’z j )l 0+ |0(1%-2j)l 0 +|DUj-(»"2) j )l 


_.N/ N .N 


< 2|D((A) r z J .)l 0+ |D(i"z j -z ;i )| 0 


(we have used the First Integral Relation for linear splines in the last inequality). 
Combining these estimates with (2.3), and noting that Dir^z= tt^Dz, we obtain (4.7). 
Finally, the convergence in (4.8) for arbitrary z e X(q) is guaranteed by (4.7), the 
dense inclusion of dom A(q) in X(q), and the uniform (in N) boundedness of 


Lemma 4.5. For z e dom A(q), 

(4,9) ||D(P N -I)u N z|| N - 0, as N-» . 

q H 


Proof: We use the triangle inequality to show 

||D(P N -I), N z|| qN < IIDIP" -iVzll^ + ||D[I N (, N z-z)]|| t|N 

+ ||D(I N z-z)|| q(| + l|D(z-A)|| qN 
< C,N|| (P N - I Vz|| „ + (c 2 + 1)||0(, N Z - z)|| „ 
+ ||D(I n z-z)|| n , 

q 

where the last inequality is due to estimates (4.5) (applicable since 
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I^A e X^(q^)) and (4.6). Thus we have 

|| D(P N - I) -rr N Z || N < C 1 N(||(P N - I)A|| n + ||A - I N A|| N ) 

+ (c 2 + 1)||(. N -I)D2|| ()N + ||D(I N z-z)|| „ 

< 2c 1 c 3 (||(^ N -I)Dz||2 n + ||D(T N z-z)|]2 )* 

* (c 2 + l)||(, N -I)Dz|| qN ♦ ||D(I N z-z)|| , 

where we have used (4.7) to argue the last inequality'. Therefore, 

the estimates in (4.1) and (4.3) may be 

used to argue || D(P^ - I)iA|| m 0 as |~| 

q N 

We now have the machinery at hand to complete our verification of the 
hypotheses of Theorem 4.1 to obtain a statement of semigroup convergence. 

Theorem 4.2. For each z e X(q) we have 

|| At;q^)iA - iA(t;q)z|| M 0 as N + «, uniformly in te[0,t], 

9r 

Proof: It remains to argue that 

||A N (q N )A - A(q)z|| « -> 0 as N 

q N 


for z e V. We see that 
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||A N (q N )A-A(.cj)z|| „ = ||P N A<q N )P N A-A(q)z|| qN 

£ l|P N A(q^(P N A - A) || ^ + ||P^(q N ) A - - N A(q)z] || qN + |I(P N -I )» N Afi)z|| qN 
£ !|A(q N )(P^A- n^z)|| qN + ||A(q^)A- A(q)z|| qN + ||(P N - I)A(q)z|| qfl 
s ei (N) + c 2 (N) + c 3 (N) . 

That e 3 (N) + 0 as N + « follows directly from (4.8). We consider e-j(N) and note that 

h(N)] 2 - | ((P N A - A) 2 (0), qi N D(P N A- A) 3 , q»D(P N A- A),)^ 

1 q 3 q^O)! (P N n M z- ir N z) 2 (0)| 2 + (q^q 2 ')|| D(P N tt N z- Tr N z)||^ 

where convergence of the second term follows from Lemma 4.5; in addition, 

| (P^z - tt^z) 2 (0) | 2 -*■ 0 as N -> « because lemmas 4.4 and 4.5 guarantee the 
H^(0,1) convergence of (P^z - ir^z) 2 . Thus e-j(N) + 0 as N -*■ «. Regarding 
e 2 (N), we use the equivalence of norms to write 

r 

t 2 (N)] 2 < v 2 ||((^ - "?)z 2 (0). (q^j --2q,)Dz 3 , (qJ.J - ,^ 2 ) Dz 2 ) T || 
2.i N N 1 2 j / n \ 1 2 . | _N N N^j i2 i n i2 

£ ]i ( | tt 2 - Tf-| | |Z 2 (0)| + |q ] Tr 3 -ir 2 qil 0 l Dz 3 L 

, i n N N N* ,2 ln , ,2, 

+ |q 2 iT 2 -Tr 3 q 2 | 0 I Dz 2 I to > , 

so that, from the convergence of -*■ 1 (i = 1,2,3), and q^ q in Q, we have 
e 2 (N) + 0asN + «. |"| 
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Just as in [7], we may use the properties of tt N , T N , and T to obtain a 

N 

statement of semigroup convergence which does not involve the tt ; similarly, 
we can obtain an analog of (4.8), i.e., strong convergence of P N to I for any 
z e X(q). Then, under the assumption that the mapping q^-*- s ( - ;q 5 ) : Ik H (0,t) 

is continuous (so that both q -* Zg(q) and q F(t;q) are continuous), we may 
obtain state variable convergence via an argument based on the "variation of 
constants" formulation of solutions. 


Corollary 4.1. For z^(*;q^) .and z ( • ;q ) solutions to (3.3) and (2.5) 
respectively, 

|| z N (t;q N ) - z ( t ;q ) |I N -*■ 0 as N -»• » 
uniformly in t e [0,t]. 

In addition, we may use the established state variable convergence to 
conclude with the desired result of this section (a slight modification of 
[5; p. 820]). 


Theorem 4.3. Assume P is compact in the P- topology (see (2.7)) and that p^ is 
a solution of (ID N ) for each N. Then {p N } contains a subsequence {p Nk } satisfying 

(i) p Nk ■*> p e P , 


N k 'vA - 

(ii) z (*;p ) + z ( * ; p ) (in an appropriate sense), and, 

(iii) p is a solution of the original parameter estimation 
problem (ID). 


Our approximation theory is not complete in that we have really only 
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considered the problem of state space approximation and have not addressed 
the problem of estimating truly variable parameters (where the infinite 
dimensional parameter space P must also be discretized). In keeping with the 
ideas of [6], [7] and [15], we shall avoid a priori parametrizations of P 
and take a more general approach that involves (linear) spline-based 


approximations for P that are independent of the level of state variable 
approximation. To this end we define P^ s I^(p) where 1^ takes elements 
p = (oiq, a-| , 6 q, 5, q-j q 4 , q^) from P and interpolates the spatially 

varying components (ct 0> a-j, Bq, B-j) to a 5- dependent mesh; that is, 

I p = (a Q , a 1( Bq, Bp C, q 3 » q 4> q g ) where, e.g., for the case of a Q , 

M M r 2M m m <\<M 

a 0 = a Q ’ = l a 0 (x k )B k (the c- dependent knots x k and linear elements B k 
k~ 0 

M 

are defined in Section 3). It is easy to see that I is continuous in the 

topology on P, so that compactness on P guarantees compactness on P^. 

We now review our convergence theory in light of these finite dimensional 

parameter spaces. Let J^(d^) = min J^(p). From the construction of P* 1 

pep" 

MM r\j M M M M M 

there exists a sequence {p ’ } in P such that p * = I p ’ ; further, using 

^i *^k N M 

compactness of P, a subsequence {p J } of {p * } may be found so that 
N i’ M k % M 

p J + p e F; using properties of I and by making additional smoothness 

JV M k * 

assumptions for P, we may argue that p d + p e P. Simple modifications 

in arguments made earlier in this section may be made to obtain corresponding 


state variable convergence as Nj,M k ». We thus obtain an analog to Theorem 4.3 

<v>N M N _M 

where a subsequence of solutions p * to the problem of minimizing J over r 

converges to a solution p to the original problem (ID). 
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5. Numerical Implementation. 

In this section we discuss some aspects of the computational algorithm for 

solving the approximate parameter estimation problem (ID N ) over the finite 

dimensional parameter space P^. To simplify the presentation we assume that 

v = 1 and that q-j = j = 1 is fixed so that an arbitrary parameter p e P has the 

form p = (3 q, e-j, 5, q 3 , q 4 , q^) with corresponding parameter in Q, 

q = (q 2 , q 3 , q 4 , q g ), q 2 (x) = Bg(x) + ^(xjH (x). In order to more easily 

describe the numerical algorithm in the case of truly variable q 2 , it will be 

more convenient to temporarily work with the parameter -7- instead of q 9 . To 

q 2 i 

this end, we note that — may be written 

q 2 


q^-(x) = Y 0 (x) + y-j(x)H ? (x) , 


where yg = 1 / gg, y-j = 1 / (Bg+B-|) - ygj in practice (when we search for a 

nonconstant q 2 ), the C[0,1] functions y Q and y-j will be estimated in place of 

Bg and ^ so we shall use the notation p = (y Q , y-j, q 3 > q 4 > q 5 ) to designate 

the unknown parameter of interest and constrain p to belong to the usual parameter 

set P. Approximations p to p, p = (yg, y-j, Z, q 3> q 4> q^) e V will be 

MM 

constructed as in Section 4, i.e., we express yg, y-j in terms of K -dependent 

linear spline elements ^ by writing y^ = £ b^ ^ for i = 0,1. The 

k-0 

approximation to then becomes, (^-) M = Yq + ^ 1* * or » 


q 2 
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C 2M w <\,ii 

l b” B"(x) , xe[0,E] 

k=0 u ’ 


(i) M (x)= { 


2M 


m ^m, 


[_ J 0 ^ b 0,k + b l,k^ B k (x ^ 


x e (§,!] , 


where S, bg ^ , and ^ (k=0, ...» 2M) are unknown and to be determined. (We 
remark that the use of Yq» y-j in place of Bg, B-j does not change our convergence 

A M 

findings due to the fact that the convergence of p to p in P yields the needed 
convergence of p M to p in P, Where now p M e (1 / Y g , -Y|/((yq + Y^Yg) > 

5, q 3 , q 4 , q 5 ) e P. ) 

N 

For given values of N and parameters, a solution z of approximating system (3.2) may be 
4N 

written z N (t;q) = £ w^(t;q)B^ where w M = col(w^) satisfies 
j =0 J J J 


(5.1) 


Q N w N (t) = K N (q)w N (t) + F N (t;q) , t e (0,t] , 


L Q n w n ( 0) - wK 


Here (F N (t ; q)) i = < F(t;q), > q and (wJJ^. = < z Q , B^ > q , for i = 0, . . . , 4N, 
while the (i,j) - elements of the (4N + 1) -square matrices and are 


given by 
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«?* 


<r n rN 
<B j’ B 1 > < 


- q 3 q 2 (o)(Bj)i(B l J , ) 1 + f (Bj) 2 (Bj )< 


0 i (B J M) 3 (B i } 3 


k n 

K ij 


= < 


A (l)Bj. B i>c 


- q 3 q 2 <0)(B^) 2 (0)CB^). 




0 q 2 


q 9 D(B^) 


i'2 




(throughout <• ,•> denotes the X(q) inner product). For a typical i and j, 
the last term in Q^. satisfies (for some 0 £ z, m <_ 2N), 

* J 


rl 

•’0 


-1- (n u ) 
qo B j 3 


( b ?>3 - 


fl 

b 


J_ rN 

q 2 


% 

B 


N 

m 


2M 


= l b 

k=0 


M 

0, k 


5 rM rN , 2 V M /.M . . M x 

B k B £ B m + ,L b 0,k + b l,k } 


k=0 


r m r n r n 
B k B * B m 


where the parameterization for approximations to l/q 2 has been used. 

'vm 

If, in contrast to the approach we take here, the spline elements B. , B , B are 

k J6 m 

defined on a uniform mesh, the appearance of? in the range of integration for 
the above integrals leads to a large amount of computational work: Such 

quadratures must be recomputed each time that 5 changes (i.e., every time the 
parameter p^ is updated in an iterative scheme to minimize over P^). An 
advantage of our formulation is that we need only compute such quadratures once. 
That is, we may use the coordinate transformation g in (3.1) to write 



f'z M N N I M N N 

where J S £ , J S^, S m need only be evaluated once, at the outset, 
then stored for recall during iteration when the coefficients 2s and 2(1- s) 
are updated. This savings in computational effort is even more substantial 
in the case of multiple discontinuities s-j, . . . , £ v , and in the case of 
unknown q-j . 

We consider now several numerical examples which illustrate the ideas presented 
thus far. In the examples that follow we return to the use of standard notation (b. will 
be used instead of ) and we assume that all parameters are known except for Bq, 6j, and S, so 
that only q£ = 3 q + ^3-j is to be determined. The reader is referred to [7] for examples with 
constant parameters in the source term s(t;qg) and in boundary conditions, and to 
[15] for multiple discontinuity examples in the context of parabolic systems. 

Indeed, both of these references provide a more complete numerical study than 
we present here: In [ 7], more realistic seismic examples are considered, 

as is the problem of surface observations (at x=0) only, while in [15], 
examples are given to illustrate that the assumption that the number of 
discontinuities is known a priori is not unnecessarily restrictive (one 
may both overestimate and underestimate the number and still get useful 
information) . 

In the examples that follow, the "true" parameter, p = (Bq> 6-j > s), is 
known and is used for comparison with our approximations, p^ = (Bq*^, , c^'). 
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Initial data for (2.1) is given, as are both velocity and pressure observations 

*\j 

(u t and q 2 u x are determined from either analytic or finite difference 
solutions for (2.1), with random noise added in some cases) at discrete time 
and spatial locations. For our numerical examples we use a "pointwise" 
fit-to-data criterion. 


(5.2) 


ftp)- i 

i » J 


C 2 r '(t.;q(p))| 

J X = X 



where C( z-j , z 2 , z 3 ) T = (0,z 2 ,z 3 ) T and y^. = (0, u t (t J .,x i ), (P 2 u x )(t J . ,x.j ) ) T ; we note 

that in using ^ (instead of the distributed criterion defined in Section 3) for our 

N 

examples, we are actually illustrating stronger convergence of z (t) to z(t) than is 
guaranteed in Section 4 (where only IR x H° x H° convergence is found). He 
thus exploit the fact that in practice IR x C x C convergence is observed, 
and only present examples that use (numerous examples with distributed 
criterion such as exist in the literature; e.g. [14] for an elliptic problem). 

We initiate the parameter estimation process by supplying an initial 
guess p^ = (Bq, B°, 5^) to IMSL's minimization routine ZXSSQ (a Levenberg- 
Marquardt algorithm) which is used in the numerical minimization of J . For 

J. l_ 

each updated value of p, the fr n approximating system (5.1) is solved using 

IMSL's DGEAR, an ordinary differential equation solver. All calculations were 

performed on the IBM 3081 D at Southern Methodist University. 

Our first example is one in which an analytic solution u is available. 

The construction of initial data and forcing function is somewhat artificial 

and serves only to guarantee an analytic solution- nevertheless, this example 

w, 

is instructive in that it is the only one in which data Cz (t) for the 
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approximate solution is compared to exact observations (Cz ( t ) ) . In later 
examples, random noise is added, or finite differences are used to solve for 
u(t) and to construct observations Cz(t). 


Example 5.1. We consider the problem of estimating the piecewise constant 
parameter q 2 with true value 



x e [0, ,.4] 
x e (.4 , 1] 


(thus = 5, (&q + $-j) = 10, % - .4). Observed data is calculated using the 
true solution 


u(t,x) 


r j (-102. 5x 2 + 96x + 48) + 12.5tx 2 , x e [0, .4] , 

.05(-31.25x 2 + 45x + 187) - ft (5x 2 - lOx + .8), x e (.4, 1] 


9 


and is available at t. = .5, 1.0, 1.5, 2.0, and x. = . 1 i , i = 0, 1, ..., 10. 

J ' 

The actual model system we use here is a nonhomogeneous version of (2.1) 

(with p = 1), i.e., 

^ ^ (E(x) f£) + f ( t , x ) , t e (0,t] , x e [0,1] , 

where boundary and initial conditions are given in (2.1). For this nonhomogeneous 
system we take f = - 9 2 u xx ’ u anc * ^2 9 1ven above, and use = 2, 

= 4, and s = 0 to construct the boundary conditions. Initial data and 
forcing function for equation (3.2) are computed using u, q 2 and f, i.e., 

Zq(x) = (u(0,0), u t (0,x), q 2 (x)u x (0,x)) and F(t,x) = (0,f(t,x), 0). 

We search for constant values of Bq and g.| , using the initial guess 



28 


^2 


'.001 , x e [0, .8] 


..001 , x e ( .8, 1] 


.N 


for N =2; we then use the converged values of q 2 to begin the iteration 
on q 2 for the next value of N. Our results are summarized in Table 5.1. 


Example 5.1 .a. We repeat Example 5.1 but add Gaussian noise to our observations 
at a level of 5% relative error. For example, if u t = average value of u t (t J .,x 1 . ) over 
all i,j, then the new data for u t is (u t ) 1 -j = u t (tj,x..) + r^., where r^. falls in the range 
[-.05u t , .05u t ] with 99% certainty. These findings are reported in Table 5.1. a. 

Example 5,1 ,b. We repeat Example 5.1 but add 10% relative noise to the 
observed data; the results may be found in Table 5.1 . b. 


Example 5.2. We present here an example that is a modification/rescaling 

(in order to include our. boundary conditions) of an example found in [ 8 ; p. 381] 

For this problem, the true value of q 2 is given by 


% 

^2 


6.25 , x e [0,^] , 


36. , x e (%,!] , 


(so (Pq, 3q + 3-|, 5) = (6.25, 36., .5)) while parameters q 2 = 1 , q^ = 6 are held 
fixed. In keeping with [8], we use Zq = (<j>(0), ip, q 2 <t> x )» w ^ ere 
<t>-(x) = exp [-160(2x - .5)^] and ip (x ) = 1600(2x- .5) exp [- 1 60 ( 2x - .5)^] , 
for initial data, and set F = 0. Observations at times t. = .05, .1, .15, .2, 

J 

and spatial locations x^ = .05i, i = 0, ..., 20, are determined by solving 

f\j 

(2.1) for u using a finite difference scheme (with q 2 , q 2 » q^, <f>, and ip 



29 




Table 5.1 : 

Example 5.1 



N 

b n 

B 0 

(3 0 + 6 1 ) N 

t" 


CP time 
(secs) 

(init) 

(0.001) 

(0.001) 

(.8000) 



2 

4.942 

8.051 

.4069 

109.05 

20. 

4 

4.959 

8.939 

.4034 

24.21 

67. 

8 

4.974 

9.382 

.4019 

5.21 

410. 

16 

4.987 

9.719 

.4009 

1.27 

1309. 

32 

4.991 

9.846 

.4005 

0.31 

1141. 

(true) 

(5.000) 

(10.000) 

(.4000) 




Table 

5.1. a: Example 

5.1. a (5% Noise Level) 


M 

a 

b n 

B 0 

<V 6 i ,N 

5 N 

^N 

u 

CP time 
(secs) 

2 

4.952 

8.064 

.4075 

109.64 

20. 

4 

4.981 

8.999 

.4042 

30.06 

65. 

8 

4.953 

8.934 

.4042 

24.87 

203. 

16 

4.959 

9.124 

.4039 

14.88 

638. 

32 

5.001 

9.877 

.4010 

7.71 

1901. 


Table 5. 

.1 .b: Example 5. 

,l.b. (10% Noise 

Level ) 


N, 

S N 

e o 

(b 0 + b 1 ) n 

c N 

j n 

CP time 
(secs) 

2 

4.986 

8.150 

.4092 

202.39 

19. 

4 

4.957 

8.001 

.4098 

194.80 

61. 

8 

4.974 

8.703 

.4071 

153.74 

185. 

16 

4.974 

8.704 

.4071 

133.75 

300. 

32 

5.072 

8.640 

.4280 

123.16 

659. 
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as above; s=0). We note that this finite difference solution for t e [0 ,.2], 
matches the graph in [8 ; Fig. 2] for t e [0 ,2.0] (our problem has been 
rescaled), and that we obtain similar wave reflection/transmission behavior at 
the interface (see Fig. 5.1). 

To estimate constant values of Sq and 8-j, we use an initial guess 
of 

f 15. , x e [0 , .4], , 



v. 15. , x e (.4,1.] , 

for N = 4 and use previous converged values to begin the N = 8, 16, 32, 40 


iterations. 

Our findings 

are given in 

Table 5.2 below. 





Table 5.2. 

Example 5.2 



N 

8 N 

B 0 

(B 0 + Bl ) N 

C N 

JN 

CP time 
(secs) 

(init) 

(15.00) 

(15.00) 

(.4000) 



4 

13.78 

16.92 

.4667 

158.09 

14. 

8 

9.78 

23.09 

.5105 

147.85 

174. 

16 

5.95 

32.10 

.4641 

60.07 

130. 

32 

6.30 

35.93 

.5000 

1.21 

260. 

40 

6.26 

35.95 

.5004 

0.44 

245. 

(true) 

(6.25) 

(36.00) 

(.5000) 





Example 5.3 Finally, we consider an example with truly spatially varying qg, 



2.5x + 1 .5 , x e [0 , .6] 

15.625X 2 - 24.375x + 9.75 ,. 


9 


x e (. 6 , 1 ] , 


(so that true c = .6). We fixed q^ = 1, q^ = 1, initial data Zq = ($(0), 

V, q 4 > x ) , where <Kx) = e and y(x) = -e , and set F = 0 ( i . e . , s = 0). Velocity and 
pressure observations were obtained by solving ( 2 . 1 ) (with the above parameter values 
and initial data) using finite difference techniques; observations 
were available at t. = .5, 1.0, 1.5, 2.0, and x. = . 05i , i =0, 1 20. 

J * 

The results of the parameter estimation process are illustrated in Figures 5.2 - 


5.4 where we compare the graphs of with those of the converged parameter (J_)^»^ 


^ V 


(we recall that - 7 - is used as the parameter in case of truly variable q ? ). Figure 5. 

q 2 

depicts the outcome of the estimation process when an initial guess of 


r 2. , x e [0,-73 

2. , x e (.7,1] 


is used, and search is made in the space p”, M = 3. 
results for M = 4 and M = 5, respectively, where the 
(M = 5) is the M = 4 (M = 5) linear interpolation of 


Figures 5.3 and 5.4 show 
initial guess for M = 4 
the converged value of 


,N,M 


/ 1 / _ \ I m 11 r u a / 11 m \ J ft 1 * n 












33 



Fiqure 5.2: Graphs of J- and (~) N,M » N - 40, M 3. 

q 2 H 2 
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